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The behavior of electrical currents in a gas discharge tube including space charge effects 
is investigated by numerical integration of the governing nonlinear partial differential 
equations. Both stationary solutions and the temporal development, under the influence 
of space charge effects, are considered. It is found that the truncation error can be greatly 
reduced by comparison with formal solutions for constant fields. The discussion is essen- 
tially restricted to the more mathematical questions. 

1. Introduction 

The bcluivior of electron and ion currents in a gas discharge tube as a function of time and 
the applied voltage has been investigated by several authors [1, 2, 3, 4]^ Most of these 
have omitted the effect of space charge, but have estimated when the effect appears. Space 
charge, however, results in a temporally growing distortion of the electrical field and, therefore, 
in a severe nonlinearity of the equations governing tlie behavior of electron and ion currents 
in the tube. A. L. Ward [5] suggested the numerical integration of the nonlinear equations 
on an electronic computer for an essentially '^one-dimensionar' tube, i.e., a tube whose elec- 
trodes are parallel plates of dimensions large compared to their separation. 

This paper describes the mathematical treatment of the basic equations necessary for the 
applicability of numerical methods. For the sake of completeness, a short derivation of the 
basic equations (sec. 2) and the formal treatment of the case of constant field (sec. 5) is given, 
although many of these considerations can be found in other papers too. The stationary case 
is treated extensively in section 3. The difference schemes used for the time-dependent case are 
discussed in section 4. The discussion of the results in section 6 is restricted to the more 
mathematical questions like the influence of truncation errors, and certain other errors occur- 
ring during the computations. A discussion of the physical significance of the results is given 
by A. L. Ward [10, 11]. 

2. Basic Equations 

We state the equations in an Eulerian coordinate system, denoting the space coordinate 
by X and the time by t. The cathode is located at x=^Q, the anode at x=d. Let n^, n^, j+, ^'_, 
v+, v^ be the density of positive ions (number of particles per unit of volume), the density of 
electrons, the current density of ions (electric charge passing through a cross section of unit 
area per unit of time), the current density of electrons, the drift velocities of ions and electrons, 
respectively. The ion current density and the ion velocity are counted positive if the ions 
are moving toward smaller x, the electron current density and the electron velocity are counted 
positive if the electrons are moving toward larger x. We denote the intensity of the electrical 
field by E, counting it positive if directed from the anode to the cathode. Let a = a{E) be 
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the number of ionizations caused by each electron per unit length of its path, q^=en-^ the 
charge density of positive ions, q^ = en^ the charge density of electrons, counted positive all 
the time, e being the elementary charge. With eo denoting the dielectric constant, the processes 
in the tube can be described by the following equations: 
(a) Continuity (Townsend^s equations) 



(2.1) 



dt 


=aj. 


" dx 


dt 


=«i- 


■^ dx 



(b) The electrostatic Maxwell's equation (Poisson's equation) 



€0 



dx 



-=g--g+. 



(2.2) 



(2.3) 



The outer circuit supplying the voltage for the tube can be described by a capacitance 
C parallel to the tube and a resistance R in series to both the tube and the capacitance. The 
outer voltage applied to this system as shown in figure 1 is called U. If V is the voltage across 
the gap of the tube and / the current to and from the tube, then the equation 



^=^(^+^f)+^ 



(2.4) 



describes the behavior of the outer circuit. The current / can be obtained from the mean 
current in the tube and the change in time of the voltage across the gap by the following con- 
siderations. Let us introduce the abbreviations q^=q~—q+ and j=j^-{-j^. Furthermore let 
S denote the area of a cross section of the current in the tube, d the distance between cathode 
and anode, and Q the charge accumulated on a unit area of the cathode. Then the law of 
conservation of charge, applied to the cathode yields 



i=s[i(o,0-f} 



Now the charge Q is connected with the voltage across the gap by 

Q=- 60^(0,0 




Figure 1. Diagram of the electric circuit. 
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and the equation 

V= f E<lx=^dE{0, t) + eo' r r ga, t)d^dx (2.5) 

Jo Jo Jo 

following from eq (2.3) by integrating twice. If we differentiate eq (2.5) w^itli respect to time 



and use r:^= — ~Lj we obtain 
ot ox 



dE(0,t) 



^Y^+ki!'^'''^''-7j^'''^- ^'-'^ 



dt 
Now the current / can be expressed as 

Introducing this into eq (2.4) leads to the following ordinary differential equation for V: 

S€n\dV 



where 



1 C^ 



is the mean current density in the gap. We write: J=qv, where v=v^v^/(v^-^v^) is an average 
velocity. Instead of eliminating dE{0,t)/dt one may eliminate dV/dt by using eq (2.6). Then, 
instead of eq (2.7), one obtains an ordinary differential equation in time for E(0,t): 

R{Cd+Seo) '^^^=U-V+^—J-R (~-+s)j{0, t) (2.8) 

at €q \ ^0 / 

where F follows from eq (2.5). 

The differential equation (2.7) or (2.8) furnishes the boundary condition corresponding 
to eq (2.3). For eqs (2.1) and (2.2) separate boundary conditions will have to be established. 
They describe the electron current at the cathode and the ion current at the anode respectively. 
The latter current is zero since there are no ions coming out of the anode: 

U(d,t) = 0. (2.9) 

The electron current at the cathode is given by 

i_(0,^)-i,+7jV(0,0+7. ('(TJ.dx (2.10) 

Jo 

where jp denotes the current density resulting from externally irradiating the cathode with 
photons. The second term describes electrons produced by ions hitting the cathode, jt being 
the probability that an incoming ion produces an outgoing electron. The last term comes 
from internally produced photons hitting the cathode. These ^^secondary^^ photons are as- 
sumed to be emitted from molecules which were excited by electron collisions. The number 
of excited molecules produced per unit of the path of a single electron is called (r=a(E), jp 
is the probability that a secondary photon produces an outgoing electron. It is assumed in 
eq (2.10) that the emission of photons occurs immediately after the excitation; however, the 
computer program contains a provision for an arbitrary time delay, simulating a delayed 
photoemission. 
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Equations (2.1), (2.2), (2.3), (2.5), (2.7), (2.9), and (2.10) together with initial distribu- 
tions of g-|_ and g;_ and the initial voltage V determine the solution of the problem completely 
for times greater than the initial time up to infinity or to a certain time limit, provided the 
velocities v^ and v^ and the ionization rates a and a are given quantities. 

Measurements show that v^, v^, and a can be approximated by the following types of functions 
of E and the pressure p in the tube: 



v,=< 



f,^{\-B,\E\/p) 



E 

V 



ii\E\<W,v 



L (Bs^^/\E\-B2f/E')E/p if \E\>W,p 
v_=fx^E/p 



a=< 



Cpexp {-D,pl\E\) 
C2P exp (-D,p/\E\) 



if \EKW2P 1 

y for molecular gases 

\i\E\>W2p J 



CsP exp {-D,^pl\E\) if \E\<JV,p ^ 



(2.11) 
(2.12) 

(2.13) 



^C,v exv {-D,^pl\E\) 



Y for atomic gases 
\i\E\>W,^ ! 



'2P J 



where jLt+, /x_, the i?„ Ci, D,, and Wi are certain constants. The quantity a is approximated by 
the same type of function as a with possibly different constants. 

3. Steady State Solutions 

The steady state has been investigated by several authors. We report here on calcula- 
tions of space charge distributions in cold cathode discharge tubes, which have been conducted 
for a number of years at the National Bureau of Standards. A. L. Ward, who suggested this 
program, reported on the results in several publications [4], [10], [11], and we shall confine 
ourselves here to stating the equations and the method to solve them. 

3.1. Differential Equations and Boundary Conditions 

We assume a state of equilibrium, which means 

dt dt 



at all times, and eq (2.1) reduces to 



^=+a(£;)i_(:r) 



(3.1) 



which is Townsend's steady state equation. In order to take account of collisions between 
electrons and metastable molecules inside the tube, we add an extra term which is proportional 

tojUx): 

'^=+a{E)j^ {x) +KE)f- (x) . (3.2) 

Here a(E) is defined by eq (2.13) and 13(E) is assumed to follow the same law, with possibly 
different parameters C„ Df, and W2. 

For eq (2.2), which governs the ions, we assume that the total current density is constant: 



J- W +j+ (^) =j = const. 
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(3.3) 



Equation (2.3) is then written as follows: 

_j(x)-j-{x) 

dE_p 



'''dx~E 



j(x)—j-(x) 



ii\E\<Wip 

(3.4) 

it\E\>W,p. 



Equations (3.2) and (3.4) are tw^o ordinary first-order differential equations lor E(x) and 
j-(x), provided that we are given the total current density function j. 

Physically accessible are the currents at the electrodes, leadino- to the followino- boundary 
conditions: 

j_(d)=j at the anode. (3.5) 

j-{0) = {j,+yu(d))/(l+y) at the cathode. (3.6) 

Here jp is a contribution to the electron current density caused by external radiation and 
7 is a secondary ionization coefficient assumed to be a constant. 

3.2. Integration of the System of Differential Equations 

We distinguisli here two cases: 
(i) The current is so small that 



dE\ E , . r» 1 1 

-7- <C-r everywhere, i.e., we assume a constant held. 

dx\ d ^ 



a(E) and P{E) are then constant also and eq (3.2) becomes a Bernoulli equation w^ith constant 
coefficients, which can be integrated explicitly. The value of E is then determined from the 
second boundary condition. 

(ii) Arbitrary large currents ^X-^)- ^^'l^^ solution lias to be found ilei-atively: 

Starting at x=d with j-{d) =j, a value of E(d) has to be assumed. (It seems most feasible 
to generate a whole family of solutions with increasing total current densities. The initial 
value E(d) is then taken to be the solution of the previous case. For sufficiently small j(x) 
assumption (i) holds and no difficulty arises in finding a starting value.) 

Equations (3.2) and (3.4) are integrated simultaneously b}^ means of a Runge-Kutta 
scheme. Iteration on E(d) is performed until j-{0) is in sufficient agreement with the pre- 
scribed boundary condition. 

As the total current densities increase, it becomes more and more difficult to find the 
proper E(d). At j_{d) ~ 10""^ amp/cm^ it was impossible to find a solution, even by the interval 
halving method. 

Fortunately this limit covers most of the experimental data as far as the basic equations 
are valid. For the results we refer to the aforementioned publication of A. L. Ward. 

4. Difference Equations 
4.1. Difference Equations With Respect to Both Time and Space. Stability Considerations 

In the most general case, it is not possible to give an exact solution of tlie system of equa- 
tions described in section 2. Therefore, one tries to find an approximate numerical solution. 
The most convenient way to obtain a ^'numerical solution^' is to introduce finite differences 
in time and space and to solve numerically the finite equations generated in this way for a 
certain set of values of the different parameters of the problem. 

The way to transform the differential equations into difference equations is largely deter- 
mined by requiring simplicity of the computational scheme and ''stability'^ of the difference 
scheme. For the sake of simplicity we ask for explicit schemes as far as possible. For the 
same reason we use difference operators of the same order as the corresponding differential 
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operators. This ensures that the special computations at the boundaries are kept to a mini- 
mum. By requiring stability we exclude certain difference schemes which would lead to 
large amplification of any small errors (as rounding errors) when the time increases, at least 
in the limit of vanishing meshwidth. We use here the concept of stability as introduced by 
Lax and Richtmyer [8, 12]. 

Denoting by A^ and Aa:; the increments in time and space respectively, the difference equa- 
tions corresponding to eqs (2.1) and (2.2) are 

for x=Ax,2Ax, . . .,MAx=d; t=tojto+At, to+2At, ... 

q^(x,t+At)-i^(x,t) ^^^^^ ^^ ._^^^ ^^^ U(x+Ax,t)-Mx,t) ^^_^^ 

for x=0, Ax,2Ax, . . ., {M—\)Ax; t=to, t^+At, to+2At, .... 

The difference quotients with respect to x have been chosen unsymmetrically for the sake of 
stability, see [8, 12]. In order to obtain stability in the sense of Richtmyer for Ai^-^0, Ax-^0, 

with — =const., we must require that 

At 

AK r^ V (4.3) 

max [v-,v^} ^ ^ 

and that both ^_ and v^ are nonnegative. Since V- is much larger than v^, this means that 
the timestep must not exceed the time an electron needs to go from one meshpoint in space 
to the next. This time is very short for many of the interesting developments in a gas discharge 
tube. Therefore, for many phenomena, it will be sufficient to assume that the electron density 
and current distribution will be in a quasi-equilibrium state, i.e., we replace the left hand side 
of eq (2.1) by zero. Since a is rather large in the interesting cases, we do not use the difference 
scheme obtained from eq (4.1) by putting the left hand side equal to zero. We rather integrate 
eq (2.1) formally: 



i_(x, 0-i-(0, t) . exp| JJa(a^', t)dx'\ 



and replace the integral in the argument of the exponential function by a finite sum, using 
the trapezoidal rule: 

{m-l 1 "\ 

T>^^[a(k^Ax,t)+aak+l)Ax,t)]Ax> (4.4) 

for m=l, 2, . . . M. 

The stability condition for the system eq (4.4), (4.2) is now 

At* 

At<~^ (4.5) 

and v+>0. Thus one can use much larger timesteps than for the original system, namely, the 
timestep must not exceed the time an ion needs to go from one meshpoint to the next. 

Equation (2.3) does not contain time derivatives. It shows that E, and therefore also a, 
are obtained by integrating q. Hence, stability is not affected, and the question how to choose 
a proper difference scheme replacing this equation can be separated from the question how to 
choose the timestep. 

Of course, the characteristic time of the outer circuit will provide another bound for the 
timestep. As long as the term containing J is unessential in eq (2.7), this characteristic time 
is apparently [R(C-\-Seo/d)]. If one uses instead of eq (2.7), the difference equation 
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R 

where 



(c+S^y^-^--^^=V{t)-V{t)-RSJ{t), (4.6) 






the stability condition for tliis scheme in the sense of Rutishauser [13] indeed is 

At<2[B{C+Seo/d)], 

as lono; as the dependence oi J on V is small enough to be disregarded. As soon as this 
dependence becomes important, we can no longer consider the outer circuit separately from 
what happens inside. Then we will have to treat the system as a whole. Even in the case 
where E and a are independent of x, i.e., when space charge effects can be neglected, this leads 
to a nonlinear problem because of the product aj_ (a depends on V). 

4.2. Difference Equation for the Electrical Field E 

The differential equation (2.3) is to be solved with the side condition (2.5) imposed on 
the integral over the unknown function E. The most reasonable way seems to be to replace 
eq (2.3) by 

E{x+Axt)-E(x, t) _ .11 



Ax 



e^'hWix+Ax,t) + q_{:x,t)] (4.7) 



for a:-0, Ax, 2Ax, . . ., {M—l)Ax. 

Using a trapezoidal rule on the left side of eq (2.5) would lead to a side condition 

V{t) = ^ /^ m{^n + \)Ax, t)+E(mAx, t))\' (4.8) 

Equations (4.7) and (4.8) form M-\-l equations for the A/+1 unknowns E(mAx, t), m = 0, 1, 

• • •' ^- 

But it turns out that these equations do not give the rigorous solutions even if q(x) is a 
linear function of x. For, take eo^q{x)=x and V=0, then, according to eqs (2.3) and (2.5), 

x'^ d^ . X? . 

E=^-^——. But, eq (4.7) yields E=—+c, where c is a constant to be determined from eq 

(4.8): 

x^ d^ (AxY 
Hence, we obtain £'==— —-^ 1-^; which is not in agreement with the rigorous solution. 

By partial integration on the right hand side of eq (2.5) one gets 

V=d.E{0,t) + eo'ij (d-x)q(x,t)dx. 
If integrated by the trapezoidal rule this formula does not give the rigorous solution for linear 

X^ rf2 ('^^\2 

Q either. For the above example we obtain E=^ + -^ — — Therefore, we have to search 

^ 2 6 6 

for a different method of replacing eqs (2.3) and (2.5) by finite equations such that the results 
are correct at least for piecewise-linear functions q. For instance, one may expect that a 
weighted mean of the two formulas discussed might eliminate the term containing (Ax)^, 
Indeed, this way is successful if the weight ratio is 2:1. The results are then correct even for 
piecewise-linear functions q, as is shown below. 

47 



Let us put €~^g[=r for simplicity. A piecewise-linear function r{x) is given by its values 
at the meshpoints rm=r{m^x) for m=0, 1, . , ., M: 

, . m^x—x , X— (m — l)Ax 
r{x)=r^-, ■ -^^+r^ ^ 

for {m — \)Lx<x<m^x {m=l,2, . . .,M). 

Our task is to find a proper approximation for the integral 

1= r rr(xyix'dx= r{d-x)r(x)dx. (4.9) 

Jo Jo Jo 

For the sake of simplicity we take d=l. The value of / follows by summation 

J*m/M 
{l—x)[r^-i(m—Mx)+7'rn{Mx—{m — l))]dx 
(m-l) 
M 

=r^., L\^ M~) 2M~6M^J+''" Lv ~mJ 2M+6M^i 
namely 

The trapezoidal rule applied to the rightmost expression in eq (4.9) gives the approximation 

J ^^ ^ f^-—\ ^ 

'' t^.MV M) 2M 

Repeated apphcation of the trapezoidal rule to the middle expression in eq (4.9) gives the 
different approximation 






^M ^0 



, 2M ' 4M' 
Hence 

l/3/„+2/37.,=[g r,„ (l-S)-?+^] h (^-^O) 

is an approximation to eq (4.9) and leads to rigorous results in the case of piecewise linear 
functions r{x) with slopes changing at x=0, /^x, 2Aa:, . . ., MAx=d=l. 

Since the solution of eq (4.7) requires the same summation as the inner sum in I^ev, we 
actually compute the finite approximation by repeated application of the trapezoidal rule and 
by adding a certain correction: 

i T 4-_ T T ^^ ^Q 

3^tr-t-3^rep-^rep 12M'* 

For general d, we have to replace r by rd^. 

The computation of E from V and the given values of g at the meshpoints can therefore 
be done by the following formulas: 

At first we compute the auxiliary quantities 

^{k-l)^x,t) form=l,2,. ..,M. 

(4.11) 



E^'imAx, t)=^ S {q(kAx, t) + q{k-l)Ax, t) for m = l, 2, . . .,M. 

^^0 k=\ 
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Then we continue with 

£'('mA/,0=-/^*(^A.^, 0+^(0,0 for m=l,2, . . .,M. (4.13) 

One may ask whether using E{{), t) instead of V{t) as the parameter describing the outer 
circuit would not simplify the computation formulas. Indeed this is true. One may compute 
£'(0, fi*oma difference approximation of eq (2.8), then compute E{mAx, t) for m=l,2, . . ., M 
from eq (4.7) successively, and finally one computes V from those values by the trapezoidal 
rule. But it turns out that the truncation error in Vis much larger in this case than if we use 
eqs (4.6), (4.11), (4.12), and (4.13). After a certain number of integration steps, F may even 
exceed t/^ which is physically impossible, when ?7= const, and C is uncharged at time ^=0. 

4.3. Method of Computation 
If we use the trapezoidal rule in eq (2.10), this boundary condition takes the form 

^ FAr "1 

i_(0,0=i.+Tri+(0,0+7.2 [Yk(mA^, ^)i_(mAx,0+^((m-l)Ax,Oi-((m-l)Aj, ^^ 

(4.14) 

The eqs (4.1), (4.2), (4.6), (4.11), (4.12), (4.13), (2.9), and (4.14) together with initial distribu- 
tions of q^ and ([^ and the initial voltage V determine the finite problem completely, provided 
the velocities v^ and v_ and the ionization rates a and a are given functions of E. Instead of 
the eq (4.1) one may use eq (4.4) if the electron density is in a quasi-equilibrium state. 

The above system of finite equations is not completely explicit. In the case where eq (4.1) 
is used, explicit formulas are achieved by taking some quantities of minor importance at an 
earlier time. Thus, in eq (4.14) the last sum has been taken at the time t—At instead of t, 
and in eqs (4.11) and (4.12), g'-(0,/) has been replaced by q^(0,t—Af). These two changes 
were sufficient to achieve explicit formulas for all quantities. 

In the electron quasi-equilibrium case the equations are even inore implicit. Therefore 
one uses an iterative scheme in order to obtain a solution starting from values of (j_ at the pre- 
vious timestep as initial approximations. The only quantity not treated iteratively is the last 
sum in eq (4.14). It is taken from the previous timestep throughout, i.e., instead of eq (4.14) 
one always uses: 

^^ At 

i_(0, t+At)=j,+y,j^{0, t+At)+y^Tj^ [a{mAx, t)j^{mAx, t) + a{{m-l)Ax, t)j^{{m-l)Ax, t)]. 

(4.15) 
4.4. Convergence of the Iteration 

It will be shown below that the proposed iteration scheme leads to convergent sequences 
for all quantities involved, if the electron density is small enough in the sense that it causes 
no distortions of the electric field comparable with the field itself. 

For the sake of simplicity we consider the limit case M-^ c» only. The quantities changing 
during the iteration are (i-{x), E(x), V-{x), a(x) void j_{x) for 0<x<d, and furthermore, ^+(0) 
and ^4.(0). All other quantities are fixed throughout the iteration. Since q-{x) is the only 
result of a previous iteration which enters the following iteration step, and since all other 
quantities depend on q-(x) continuously, it is sufficient to prove that the sequence of functions 
q-{x) converges. 

According to [7], one has to consider the change dq^'j'^ix) of the result of one iteration 
step, caused by a certain change 6g!!^ (x) of the initial approximation. When measured by a 
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certain norm, the ratio of the changes has to be smaller than unity in a certain neighborhood 
of the true solution q^(x). 

Tlie following formulas describe the connection between 8q^^'' (x) and dq''}^ (x). They are 
derived from the continuous analogs of the eqs (4.4), (4.11), (4.12), (4.13), from (4.15) and 
from the relations between v^, ^+, a and E, and those between g±, v^^ and j^. 






i-(o) 

and r is the last term on the right-hand side of eq (2.10). 

8v-(x)_ 1 dv 



V- (x) v- (IE 



8E(x) 



da(x)=^dE(x) 



8E{x)=eo' \ K{x,u)bq^^{u)d^ 
Jo 

K{x,u)-- 



Jo 
where 

u/d' if u<^x 

ujd— 1 if u>x. 



We restrict ourselves to the case for which not only v_ but also z;+ is proportional to E. Then 

l_dvj^__ 1 dv-_ 1 
v^dE~vldE~'E 
Furthermore, we write 

da a (/(In a) 

dE~Ed{\nE) 

since the latter differential quotient varies more slowly in the interesting range of E. We 
introduce the norm 

\\j\\=^l\m\dx 

for any function /(x). Then, the following estimates can easily be derived: 



|5£:(x)|^e„-'max||l-^|.||52- 

€o> r \<p^{:x)\dx 

\\5q^-\\^ K '^'. ,„, ,, ■||3gi"'|| 
min \E{x)\ " ^ " 



X 

where 

min|£(x)| 3 

^=1+ i?(0) +4"'^^?^ ^^(ta^ 



The condition for converaence is that 



€o-' ( \q'^^{x)\dx 
^ Jo 



K -^^ , - ^. ., ^Q<1 (4.16) 

mm |A(ic)| 
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wliieli can be interpreted the following way: The field distortions produced by the electrons 
and aini^lified by the numerical factor ii must not exceed the minimum field strength in the gap. 
From the rate of convergence given by eq (4.16) and the change of q-{x) between suc- 
cessive iterations, tlie corrections introduced by a further iterative step can be estimated. 
According to these considerations, the iteration may be stopped after a prescribed accuracy 
has been reached. A criterion of tliis kind has been used in the code. 

5. Formal Solutions for Constant Field 
5.1. General Considerations and Formulas 

The nonlinearities contained in the equations of section 2 and in the corresponding differ- 
ence equations disappear as soon as the quantities a, a, Vj^ v- can be regarded as independent 
of the solution. If, moreover, these quantities are constant in space and tune and if the voltage 
V across the gap is constant, the equations of section 2 and the corresponding difference 
equations can be solved explicitly. Then, of course, there is no room for eq (2.3), i.e., this 
treatment disregards space charge, andeq (2.7) cannot be taken into account, i.e., the reactions 
in the outer circuit are disregarded, or, in other words, the external resistance R is so small 
that it can be neglected by putting 7?=0, which leads to U=V. Work in this direction has 
been done by several authors [1, 2, 3, 9] as far as the differential equations are concerned. These 
papers discuss what happens in the tube. Here this special case will be considered again, but 
for a different purpose. We shall discuss the difference equations along with the differential 
equations in order to get insight into the effects of truncation errors. We hope that, to a 
certain extent, these effects carry over to the more general case, and therefore will allow 
us to correct the results obtained with finite steplengths A^ and Ax so that we obtain closer 
approximations to the case of infinitesimal steplengths. 

We mainly deal with the difference equations. The results for the differential equations 
will be obtained by letting A^, Aa^->0. We restate the equations of section 4 for our special 
case putting ^o = 0: 

j^(x,t+At)-j-(x,t) ^., s j^(x,t)-j-(x-Ax,t) . . 

V-At J \ ^ J ^^ K • <) 

for x=Ax, 2Aa:, . . . MAx=d] t=0, At, 2At, . . . 

In the electron quasi-equilibrium case, we shall use instead of this: 

Mx,t)=j-(0,t) .e«^ (5.1b) 

for the same x and t as above. The ion currents behave according to 

j+(x,t+At)-j+(x,t) _. , . j+(x+Ax,t)-j^(x,t ) 

^^ --j_[x,t)^ ^^ [5.Z) 

for x=0, Ax, , . . (M-l)Ax; t=0, At, 2At, . . . 
There are two boundary conditions, namely 

i+(rf, t) = for ^-0, A^, 2A^, . . . (5.3) 

and 

j.[0, t)=j,+y,MO, 0+7. 2 ^ a-[j.(mAx, t-At) +j_((m-l)Ax, t-At)] 

m=l ^ 

for t = At, 2 At, . . . (5.4) 

In order to make the solution unique, we have to in troduce initial conditions, for instance by 
prescribing j-{x,0) smd j+{x, 0) for x=0, Ax, , , , MAx. But we will be mainly concerned with 
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solutions which are proportional to e^^ with a suitable X. For this type of solution, initial 
conditions do not have to be given. 

The above equations are linear and homogeneous with the single exception of (5.4) which 
contains the inhomogeneous term j^j. Since jjy is time independent and since time does not 
appear explicitly in any other term, there will exist, at least in general, a time independent solution. 
More general solutions will be obtained by superposition of any solutions of the corresponding 
homogeneous system with j'^=0. The homogeneous system allows for time separation in the 
form of a factor 6^^ Then the stationary case results from putting X=0 and from slight changes 
due to the inhomogeneous term. We introduce jx^ {x) and j^;" {x) by 

i_(x, t)=j^{x) -e^' 

(5.5) 
j+{x, t)=j^(x) -e^' 

Then, eq (5.1a) can be solved by 

ixW--ix(0)•6«*^^ (5.6) 

where a* is the solution of the transcendental equation 



e^^' — l l-6-«*^^ 



v^'At Ax 



(5.7) 



In the electron quasi-equilibrium case, i.e., if we use eq (5.1b), all results will be correct 
in the sequel if we replace a* by a, unless we distinguish explicitly between the two cases. 
Equation (5.2) now reduces to a single inhomogeneous equation ior j^(x)y the inhomogeneous 
term being ajx{x). The solution satisfying the boundary condition (5.3) is easily found to be 

j^(x)=Ae^'[e^<^*-^^^-e^-'-^^''], (5.8) 

where 

h(0) 



aAx V-^aAt 

and ^ is determined by the equation 

v+At Ax 



(5.9) 



The remaining unknowns are j^^ (0) and X. In the homogeneous case the quantity j^;; (0) remains 
free, whereas X follows from eq (5.4) with jp = 0, which transforms into the transcendental 
equation 

aAa; v^aAt 

5.2. Stationary Solution 

For the stationary solution of the inhomogeneous eq (5.4), we put X=0. Then (5.9) 
leads to ^=0, and (5.4) shows that 

•-^m • //i Q^Ax(6-*"-l) / , (7H-6«^\^ . ^ - 

.70 {o)=j,l^i — ^a*Ax_i \yi+yp- — 2 — ) r (5-11) 
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If we use (5.1a), we may simplify our equations further. From eq (5.7) we obtain a* 
explicitly: 

«*=;r-liir-^- (5.12) 

Ax 1— aAx 
Equation (5.8) reduces to 

J^(x)=Jom • (l-cxAx){e^''-e-'^), (5.13) 

and eq (5.11) reduces to 



io-(o) 



==i.//l-[^"*'-l][7.(l-c.Ax)+7.^(l~|^Ax)]|^ (5.14) 



In the electron quasi-equilibrium case, we use eq (5.1b) and obtain 

h^ (x) =i„- (0) • ^^ ie""- e"^) (5.15) 



and 

io-(o)=i 



./{'-^ <'■'-» (-+-;'^')}- 



By letting Aa:-^0 we obtain the corresponding relations for the differential equations; 
namely a'^ = a and, furthermore, 

ioMa;)=io-(0)-(«''^-0, (.'5.17) 

io-(0)=i./{l-(<^"^-l) (t.+7„^)}- (5.18) 

When the respective denominators in eqs (5.14), (5.16), (5.18), become zero, no stationary 
solution is possible. The voltage V for which tliis occurs is called the breakdown voltage. 
If we assume that aja is a given constant independent of V (it is often assumed that (j = a) then 
the three respective denominators are equal to 1 for a = and decrease monotonically in the 
interval ()<a<i^ . Therefore, there exists one and only one positive real root a, corresponding 
to the breakdown voltage. If the voltage V exceeds this value, the tlieoretical stationary 
solution shows negative values of y~, which is physically impossible: No stable solutions exist 
beyond the breakdown voltage. 

The physically interesting things happen if V is near the breakdown voltage. Therefore 
we have to look into the dependence of the breakdown voltage on the steplength Ax in order 
to establish a base for comparison of numerically obtained results and the true theoretical solu- 
tions, or the experimental measurements in this case. 

In the infinitesimal case the value a = a^^ corresponding to breakdown can be expressed 
explicitly from (5.18), if a I a is given: 

For the difference schemes, the corresponding transfornnition will lead to a form of the equation 
suitable for iteration. From eq (5.16), i.e., for the electron quasi-equilibrium case, we get 

In case we use eq (5.1a) instead, it follows from eqs (5.14) or (5.11) that 

„v=i.„[i+»p«.,)/(„+..(^X_i±^5V^)]. 
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We shall see that these formulas are special cases of more general formulas for the homogeneous 
system. 

5.3. Special Time-Dependent Solutions 

We turn now to a discussion of the homogeneous system, especially the transcendental 
eq (5.10). Let us first consider the infinitesimal analog of this equation, obtained by letting 
Ax^O and At-^0, Equation (5.9) becomes then 

/5-XK, (5.19) 

and eq (5.7) transforms into 

a^=:,a-\/v^- (5.20) 

If we define z; by 1/^; = 1/v^+l/v^, then eq (5.10) becomes 

a — K/V a — K/V- 

in accordance with [9]. For the breakdown voltage, X=0 is a solution. 

For voltages near the breakdown voltage, there must be a solution X near zero for continuity 
reasons. For this special root, the equation may be transformed into an iterative scheme similar 
to the one for the computation of breakdown voltage itself. We show this for the case where 
a solution is wanted for X as a function of a. 
We put 

a—X/v=a(l — u) 
and rewrite eq (5.21): 

According to the assumptions made above, only the first factor on the right changes rapidly 
with a. Therefore we solve the equation for this factor, obtaining 

This form is suitable for iteration since the right hand side is only slowly varying with a, par- 
ticularly, if we use u as an independent parameter. We describe the function \{a) by the 
parametric representation \=X(u), a=a(u). From a first approximation of V and awe com- 
pute v^, V-, V, a, and 

\ = avu. (5.23) 

Since all variables on the right hand side of eq (5.22) are known, a new approximation for a 
can be computed from eq (5.22). Since the right hand side is slowly varying with a for suffi- 
ciently small Uy we can hope for fast convergence of the iteration. 

A similar but slightly more complicated approach is possible for solving eq (5.10), if eq 
(5.1a) is used for the computation of electron current. Here, in addition, the relation between 
a* and a has to be established iteratively. Let us assume that we know a first approximation 
for a* for a given value of the independent parameter 

"^^-^ — rr" 

avAt 

Then a can be computed from (5.7) if we assume in addition that v/v^ is at least approximately 
known: 
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Ax 



'('-")■ 



(5.24) 



If necessary, from the a thus ol)taino(l one may compute a new value of viv- and repeat the 
computation of a iteratively. 

From u and a one ma^^ compute 



and 



X=— hi (l+avuAt) 






(5.25) 



(5.26) 



Equation (5.10) may be written as 



l=[e^a*-m-i]. 



ga*Aa: ( 1_^ — )-t^ — 6«*^^ (\—U — ^ 

V V.J v+ V V.J 

and is solved for the first factor: 



2 * i_e-(«*-^>'^ 



1 



(1-^/a*)^ 



In 



1 + / 



Tz 



7;,(a/a;)e-^^^ 1 + 6"*^^ 6' 



,)3d z,-(a* -/3)d\ -1 



.a*A. A _^ ^)-^ ^ ^"*^' (l-^ ^) 



]^ g-(a*-/3)d 



(5.27) 

The computation of the breakdown voltage thus appears as a special case of these formulas, 
namely as the case u=\=0. 

In the equilibrium case the formula corresponding to eq (5.27) is 



1 



"(l-^/a)r/ 



In 



1+/ 



7e 



yj,{(jla)e^ 



{e<^^^-l)l{aAx)-u 



V ' (6«^^-1)/(q;Ax) 



2 l — (i-ia-fi)d I 



(5.28) 



This equation is used together with eqs (5.25) and (5.26) for solving eq (5.10) iteratively. 

Because of the similarities between the eqs (5.22), (5.27), and (5.28) one can treat all 
three of these equations quite simply in a single computer program. 

5.4 General Time-Dependent Solution for Constant Field 

So far we have considered only a special root of eq (5.10). However, one can see that a 
transcendental equation of the type of eq (5.10) and its limit case eq (5.21) has more roots X 
in general. For the case of finite A^ and Ax, the equation is rational in e^^K Therefore, there 
is a finite number of roots e^^^, each of which corresponds to a infinite set of roots X of the 
form \=\Q-{-2Tik/At (k=0, ±1, ±2, . . .). But all these roots describe the same function 
on the grid. Equation (5.21), however, has an infinity of roots X, whose asymptotic distribution 
is shown in [6]. It turns out that at most a finite number of them can have a real part larger 
than the real root discussed before. The root with the largest real part will become dominant 
in any solution as time goes on. Therefore, one is mainly interested in the root with the 
largest real part. If we knew that this root was the real root discussed in section 5.3, we 
could confine our considerations to this root essentially. Unfortunately, a proof is not avail- 
able at the time being. Therefore, it remains an open question whether or not the asymptotic 
behavior of the solution for larger t can be described by the formulas given in section 5.3. 
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For a complete formal solution of an arbitrary initial value problem, of course, one has to 
consider all of the solutions of eqs (5.10) and (5.21) and one has to develop the initial distri- 
bution into a series of the corresponding functions after substracting the stationary solution 
of the nonhomogeneous equation. No attempt has been made to go further in this direction, 
but see [1] for some results of this kind. 

6. Results 

The physical significance of the results obtained by the machine computations is discussed 
in [5, 11]. Therefore, we restrict our discussion to the more mathematical questions. The 
main question is how large the truncation error is, i.e., the error introduced by using finite 
differences instead of derivatives. No attempts to establish rigorous error bounds have been 
made. Instead of this, experiments with different step lengths have been carried out for the 
following set of parameters: 

jp=10~'^^ amp/cm^, p = 722 torr, d=l cm, aS'=1 cm^, 

'C=10-'' amp sec/volt, R=10^ volt/amp, U=V=25.6 kv for /<0. 

U for ^>0, ji, 7p and Ax are different for different curves. The time step At is given by the 
formula A^= 0.8- Ax/max {?;+} in the electron quasi-equilibrium case, A^=0. 8* Ax/max {vJ] in the 

X X 

general case, unless something else is stated explicitly. The constants used to describe the 
functions a. and a are Ci = C2=^ cm~^ torr~\ Di=D2 = 24:7 volt cm~^ torr~\ the mobilities of 
electrons and ions are constant: jli_=4X10^ torr cmV(sec volt), ijl^ = 2X10^ torr cm7(sec volt). 

The calculations in section 5 show that in the most interesting area near the breakdown 
voltage the solution is very sensitive to changes in the voltage. On the other hand, a finite 
stepwidth of reasonable size, e.g., with 20 subintervals, introduces a change of the breakdown 
voltage of notable magnitude. Therefore, the differences between runs with different step- 
widths are mainly due to the influence of the stepwidth on the breakdown voltage. It seems 
to be feasible to eliminate this influence by relating the applied voltages to the breakdown 
voltage, as computed according to section 5.3, i.e., neglecting space charge, for the stepwidth 
used in each case. This method turned out to give very satisfactory results for the electron 
quasi-equilibrium case, as can be seen from figures 2 and 3. 

The reason why even the relatively large steplength Ax=0.05 cm (a:Ax«0.5 to 0.6) gives 
a good approximation, can be seen from the figures 4, 5, and 6, where the coefficient of temporal 
growth X, as computed from the formulas of section 5.3, i.e., without space charge effects, is 
plotted versus the voltage V across the gap and the overvoltage V— V^j,. Figure 4 shows that 
the curves for the difference equation (quasi-equilibrium case, Ax =0.05 cm) and for the differ- 
ential equation go almost parallel over a long range. Therefore after relating the voltages 
to the breakdown voltage, the curves almost coincide as can be seen from figure 5. This 
explains the good results obtained with that approximation, at least as long as space charge 
effects have small influence. 

With the same value of ji+yp, but a portion of 10 and 20 percent jp, the temporal devel- 
opment goes faster than for 7^ = 0, as one could expect from physical considerations. 

According to flgure 6, the difference scheme should give an approximation almost as good 
as for yp=0, even a better one in the 10 percent case, where the curves for the differential 
equation and the difference solution nearly coincide. 

The influence of the time lag in the term describing the production of electrons by secondary 
photons (last term in eq (5.10)) has been studied by introducing an artificial factor exp(— XA^) 
with jp into eq (5.22) for the continuous case. A^ was assigned a fixed value approximately 
equal to the ones used in the computation with Ax=0.025, namely A^=6. 25X10"^ sec. The 
deviations due to that factor exp(— XA^) can be seen from table 1 below. 
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Figure 2A. Mean current density J across the gap versus 
time t for yp = 0, yi= 1.5 10"^ 

Ax [cm] U[kv] C7-Fbr[kv] 

(a) 0.05 28.870 0.233 

(b) .025 28.870 .295 

(c) .025 28.807 .233 



Figure 2B. Mean current density J across the gap 
versus time t for 7p=1.5 10~^ 7i = 1.35 10~^. 



x[cm] U[kv] C7-Ubr[kv] 

(d) 0.05 28.854 0.233 

(e) .025 28.854 .287 

(f) .025 28.800 .233 



We conclude that the influence of tliat time lag is not very important, at least for the 
overvoltages and the small rates 7^ of production of secondary electrons considered here. 
We see that the time lag slows down the speed of development by a few percent at most, even 
in the worst case. 

The electron quasi-equihbrium assumption dq„/dt=0 in (2.1) is equivalent to letting 
?;_-^oo in this equation. Therefore, the influence of that assumption in the differential equa- 
tion may be studied by replacing v^ by 00 and v by v+ in the eqs (5.22) and (5.23). If yp=0, 
the result is obvious. For fixed voltage, i.e., fixed a, the introduction of that assumption 
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Figure 5. Coefficient X of temporal growth for coiistant field 
versus overvoltage V — Vbr for same parameters as figure 4- 
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Figure 3A. Voltage V across the gap versus time t for param- 
eters of figure 2A. 

Figure 3B Voltage V across the gap versus time t for param^ 
eters of figure 2B. 
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Figure 4. Coefficient X of temporal growth for constant field 
versus voltage V across the gap, for Ax = 0.05 C7n. 

(a) 7p=0, 7i=i.5 10-5 

(b) 7p= 1.5 10-6, 7,= 1.35 10-5 (contmuous case only) 
for electron quasi-equilibrium ease (dashed) 

general case (dotted) 
contmuous case (solid). 
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Figure 6. Coefficient X of temporal growth for constant field 
versus overvoltage V — Ybr for Ax = 0.05 cm. 

(a) 75=0, 7t=1.5 10-5 

(b) 7p=1.5 10-8, 7i=i.35 10-3 
(continuous and quasi- 
equilibrium case coin- 
cide) 

(c) tp=3.0 10-6, 7i=1.20 10-5 

for electron quasi-equilibrium case (dashed) 
continuous case (solid). 
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Table 1. Influence of the time lag Ai for secondary photons on the coefficient X of temporal growth 



F-Fbr 


T" -0.1 

Tj+Tp 


T: -0.2 
7j+7p 


A< = 


A«=6.25 10-7 


A^ = 


i = 6.25 10-7 


VoU 
-400 
-300 
-200 
-100 

+100 
+200 
+300 
+400 


-0. 8643 
-. 6595 
-. 4479 
-. 2285 


+0. 2396 
+. 4930 
+. 7605 

+1. 0610 


-0. 8624 

-.6576 

-.4463 

-.2274 



+0. 2377 
+. 4878 
+.7538 

+1.0410 


-0. 9291 

-. 7147 

-. 4901 

-. 2532 



+0. 2746 
+. 5802 
+. 9317 

+1. 3650 


-0. 9248 
-. 7105 
-. 4864 
-. 2506 


+0. 2696 
+. 5651 
+. 8978 

+1. 2880 



amounts to no change in u, and therefore to a relative increase of X by about v^lv_. Hence 
the temporal development is speeded up by about 0.5 percent, which is small enough for the 
accuracy required in this problem. 

A detailed analysis of the effect of that assumption in the case t^f^O, which is elementary, 
but too lengthy to be reproduced liere, shows that this remains true as long as the relative 
influence of y^ as compared to 7^ in eq (5.22) is small, i.e., roughly speaking, aslongas7p-6^'^^''+<^7i. 
This is in agreement with physical considerations, since, if the influence of secondary 
photons is dominant, the development is dominated by a process with feedback, all of whose 
components go infinitely fast. Quantitatively, as long as ti^I, v+<^v^, ad^l, the change 
8\ of X effected by the equilibrium assumption is, with good approximation. 



X v^ \ 



ad— I yj,{ala){\—u)e^^"'^ 



ad-ll{l-u) 



The full difference equations, with eq (5.1a) included, have not been used^, for the follow- 
ing reasons. For the same Ax and the same time interval to be covered, the stability conditions 
require 200 times as many time steps, if the electron-quasi-equilibrium condition is dropped. 
Furthermore, figures 3 and 4 show that, for Ax=0.05 cm, the deviation from the limit case 
^x=Q is considerable, and cannot be diminished satisfactorily by a simple change in the outer 
voltage leading to the same overvoltage. It is estimated that a decrease of A^:; to 0.01 cm is 
at least necessary in order to obtain results as good as the ones in the quasi .equilibrium case. 
Two trial runs were made, both leading to breakdown after a few ion transit times. The 
data for these runs were: 

(a) Ax=0.05, ^7=31.5 kv, Fbr = 28.64 kv, 7z=1.5 10"^ 7;,=0 
infinite current after 2 ion transit times. 

(b) Aa:=0.5, ?7- 28.87 kv, V^,=2^.bl kv, 7z=0, 7^=1-5 10"^ 
infinite current after 3 ion transit times. 



All the above considerations and the computations done so far are restricted to a domain 
where the temporal growth in one (natural) time-step (as indicated by the stability condition) 
is not too large. Near the actual breakdown, therefore, these arguments may not apply, 
and an investigation of the nature and quantitative features of the breakdown singularity 
is not given in this paper. Since the time required for the space change to become significant 
is very large compared to the remaining time until completion of the breakdown, and since 
the assumption of one-dimensionality is doubted for high current densities, there was no 
point in trying to describe the details of the breakdown more precisely. 

The features discussed in this section have been tested only for a few sets of parameters. 
Up to this time it has not been proved that the conclusions are general. Since the number of 



2 Several successful computer runs have been made since by Dr. A. L. Ward, using this option of the program. 
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parameters is quite large, it has not been attempted to explore the limits of the region where 
the conclusions are valid. But the methods of computing and of guessing the effect of trunca- 
tion errors as described here should be tried in other cases. 



The authors express their thanks to Dr. A. L. Ward of the Harry Diamond Laboratories, 
Department of the Army, for suggesting the problem, and to J. P. Menard and A. E. Beam, 
who did the computer programming. 
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